Appendix B: Agricultural Insurance Loss and Relationships to Climate across the Inland Pacific Northwest Region of the United States
Erich Seamon, Paul E. Gessler, John T. Abatzoglou, Philip Mote, Stephen Lee
September 2019
Appendix B Summary Statement: Agricultural Insurance Loss and Principle Components Analysis
Appendix B documents the exploratory data analysis process, examining the relationships of agricultural commodity loss, at a county level, from 1989-2015, across the three state region of the Pacific Northwest (Washington, Idaho, and Oregon), and then focus in on the 24 county region of the Inland Pacific Northwest (IPNW).
Here we perform a variety of Principle Components Analyses (PCA) to better understand the combined effects of differing damage causes, commodities, counties, and years on overall loss.
Step 1. Data Preparation. Here we load our agricultural insurance loss data and prepare it for PCA Analysis.
Step 2. Data Preparation. Original Insurance Loss Dataset - Pacific Northwest
Step 3. Data Preparation. Aggregated Insurance Loss Dataset - Pacific Northwest
Step 4. 24-county inland Pacific Northwest (iPNW) Study Area.
Step 5. PCA: PNW insurance loss, by COUNTY with COMMODITY loadings: 2001 to 2015
Step 6. PCA: PNW insurance loss, by YEAR with DAMAGE CAUSE loadings: 2001 to 2015
Step 7. PCA: PNW insurance loss, by MONTH with DAMAGE CAUSE loadings: 2001 to 2015
Step 8. PCA: PNW WHEAT insurance loss, by YEAR with DAMAGE CAUSE loadings: 2001 to 2015
Step 9. PCA: PNW WHEAT insurance loss, by COUNTY with DAMAGE CAUSE loadings: 2001 to 2015
Step 10. PCA: PNW APPLES insurance loss, by COUNTY with DAMAGE CAUSE loadings: 2001 to 2015
Step 11. PCA: PNW BARLEY insurance loss, by COUNTY with DAMAGE CAUSE loadings: 2001 to 2015
Step 12. PCA: IPNW insurance loss by COUNTY with DAMAGE CAUSE loadings: 2001 to 2015
Step 13. PCA: IPNW insurance loss, by COUNTY with DAMAGE CAUSE loadings: 2001 to 2015
Step 14. PCA: IPNW insurance loss, by YEAR with DAMAGE CAUSE loadings: 2001 to 2015
Step 15. PCA: IPNW insurance loss, by COUNTY with DAMAGE CAUSE loadings: 2015
Step 16. PCA: IPNW WHEAT loss by YEAR with DAMAGE CAUSE loadings: Temporal Mapping of PC1 and PC2
Step 17. PCA: IPNW WHEAT loss by COUNTY with DAMAGE CAUSE loadings: Spatial Mapping of PC1 and PC2
Step 1: Data Preparation
The first step was to load Pacific Northwest agricultural insurance commodity loss data, aggregate by county, damage cause and commodity, remove zeros, transform the loss values ($) into cube root and logrithmic outputs, and scale and center the data.
Original data file for all insurance claims for the Pacific Northwest, can be found here:
https://github.com/erichseamon/paper-PHD1/blob/master/data/RMA_PNW_2001_2015.csv
Aggregated data file, which summarizes data by year, county, commodity, and damage cause, for just the Pacific Northwest, can be found here:
https://github.com/erichseamon/paper-PHD1/blob/master/data/RMA_damage_PNW_2001_2015.csv
All code that generates the following graphs and tables, can be found here:
https://github.com/erichseamon/paper-PHD1/
library(tab)
library(car)
library(RCurl)
library(lme4)
library(ez)
library(lattice)
library(ggplot2)
library(coefplot2)
library(broom)
library(htmlTable)
library(gridExtra)
library(kableExtra)
library(ggfortify)
#options(scipen=999)
#load original data as txt file
rma1988 <- read.csv("/dmine/data/USDA/crop_indemnity_txt/1988.txt", sep = "|", header = FALSE, strip.white = TRUE)
rma1989 <- read.csv("/dmine/data/USDA/crop_indemnity_txt/1989.txt", sep = "|", header = FALSE, strip.white = TRUE)
rma1990 <- read.csv("/dmine/data/USDA/crop_indemnity_txt/1990.txt", sep = "|", header = FALSE, strip.white = TRUE)
rma1991 <- read.csv("/dmine/data/USDA/crop_indemnity_txt/1991.txt", sep = "|", header = FALSE, strip.white = TRUE)
rma1992 <- read.csv("/dmine/data/USDA/crop_indemnity_txt/1992.txt", sep = "|", header = FALSE, strip.white = TRUE)
rma1993 <- read.csv("/dmine/data/USDA/crop_indemnity_txt/1993.txt", sep = "|", header = FALSE, strip.white = TRUE)
rma1994 <- read.csv("/dmine/data/USDA/crop_indemnity_txt/1994.txt", sep = "|", header = FALSE, strip.white = TRUE)
rma1995 <- read.csv("/dmine/data/USDA/crop_indemnity_txt/1995.txt", sep = "|", header = FALSE, strip.white = TRUE)
rma1996 <- read.csv("/dmine/data/USDA/crop_indemnity_txt/1996.txt", sep = "|", header = FALSE, strip.white = TRUE)
rma1997 <- read.csv("/dmine/data/USDA/crop_indemnity_txt/1997.txt", sep = "|", header = FALSE, strip.white = TRUE)
rma1998 <- read.csv("/dmine/data/USDA/crop_indemnity_txt/1998.txt", sep = "|", header = FALSE, strip.white = TRUE)
rma1999 <- read.csv("/dmine/data/USDA/crop_indemnity_txt/1999.txt", sep = "|", header = FALSE, strip.white = TRUE)
rma2000 <- read.csv("/dmine/data/USDA/crop_indemnity_txt/2000.txt", sep = "|", header = FALSE, strip.white = TRUE)
rma2001 <- read.csv("/dmine/data/USDA/crop_indemnity_txt/2001.txt", sep = "|", header = FALSE, strip.white = TRUE)
rma2002 <- read.csv("/dmine/data/USDA/crop_indemnity_txt/2002.txt", sep = "|", header = FALSE, strip.white = TRUE)
rma2003 <- read.csv("/dmine/data/USDA/crop_indemnity_txt/2003.txt", sep = "|", header = FALSE, strip.white = TRUE)
rma2004 <- read.csv("/dmine/data/USDA/crop_indemnity_txt/2004.txt", sep = "|", header = FALSE, strip.white = TRUE)
rma2005 <- read.csv("/dmine/data/USDA/crop_indemnity_txt/2005.txt", sep = "|", header = FALSE, strip.white = TRUE)
rma2006 <- read.csv("/dmine/data/USDA/crop_indemnity_txt/2006.txt", sep = "|", header = FALSE, strip.white = TRUE)
rma2007 <- read.csv("/dmine/data/USDA/crop_indemnity_txt/2007.txt", sep = "|", header = FALSE, strip.white = TRUE)
rma2008 <- read.csv("/dmine/data/USDA/crop_indemnity_txt/2008.txt", sep = "|", header = FALSE, strip.white = TRUE)
rma2009 <- read.csv("/dmine/data/USDA/crop_indemnity_txt/2009.txt", sep = "|", header = FALSE, strip.white = TRUE)
rma2010 <- read.csv("/dmine/data/USDA/crop_indemnity_txt/2010.txt", sep = "|", header = FALSE, strip.white = TRUE)
rma2011 <- read.csv("/dmine/data/USDA/crop_indemnity_txt/2011.txt", sep = "|", header = FALSE, strip.white = TRUE)
rma2012 <- read.csv("/dmine/data/USDA/crop_indemnity_txt/2012.txt", sep = "|", header = FALSE, strip.white = TRUE)
rma2013 <- read.csv("/dmine/data/USDA/crop_indemnity_txt/2013.txt", sep = "|", header = FALSE, strip.white = TRUE)
rma2014 <- read.csv("/dmine/data/USDA/crop_indemnity_txt/2014.txt", sep = "|", header = FALSE, strip.white = TRUE)
rma2015 <- read.csv("/dmine/data/USDA/crop_indemnity_txt/2015.txt", sep = "|", header = FALSE, strip.white = TRUE)
#load original data as from github as a joint file - 1989 to 2001 and 2001 to 2015
#RMA_2015 <- read.csv("https://raw.githubusercontent.com/erichseamon/paper-PHD1/master/data/RMA_csv/US_RMA_2001_2015.csv", sep = ",", header = FALSE, strip.white = TRUE)
#RMA_1989 <- read.csv("https://raw.githubusercontent.com/erichseamon/paper-PHD1/master/data/RMA_csv/US_RMA_1989_2000.csv", sep = ",", header = FALSE, strip.white = TRUE)
RMA_1989 <- rbind(rma1989, rma1990, rma1991, rma1992, rma1993, rma1994, rma1995, rma1996, rma1997, rma1998, rma1999, rma2000)
RMA_1989 <- data.frame(RMA_1989[,1], RMA_1989[,3], RMA_1989[,5], RMA_1989[,7], RMA_1989[,12], RMA_1989[,13], RMA_1989[,14], RMA_1989[,15])
colnames(RMA_1989) <- c("year", "state", "county", "commodity", "damagecause", "monthcode", "month", "loss")
#write.csv(RMA_1989, file = "/dmine/code/git/paper-PHD1/data/RMA_csv/US_RMA_1989_2000.csv")
#save(RMA_1989,file="/dmine/code/git/paper-PHD1/data/RMA_Rda/US_RMA_1989_2000.Rda")
RMA <- rbind(rma2001, rma2002, rma2003, rma2004, rma2005, rma2006, rma2007, rma2008, rma2009, rma2010, rma2011, rma2012, rma2013, rma2014, rma2015)
RMA <- data.frame(RMA[,1], RMA[,3], RMA[,5], RMA[,7], RMA[,12], RMA[,13], RMA[,14], RMA[,15], RMA[,16])
colnames(RMA) <- c("year", "state", "county", "commodity", "damagecause", "monthcode", "month", "acres", "loss")
#write.csv(RMA, file = "/dmine/code/git/paper-PHD1/data/RMA_csv/US_RMA_2001_2015.csv")
#save(RMA,file="/dmine/code/git/paper-PHD1/data/RMA_Rda/US_RMA_2001_2015.Rda")
RMA_PNW <- subset(RMA, state == "WA" | state == "OR" | state == "ID" )
#save(RMA_PNW,file="/dmine/code/git/paper-PHD1/data/RMA_Rda/PNW_RMA_2001_2015.Rda")
#write.csv(RMA_PNW, file = "/dmine/code/git/paper-PHD1/data/RMA_csv/RMA_PNW_2001_2015.csv")
RMA_PNW_1989 <- subset(RMA_1989, state == "WA" | state == "OR" | state == "ID" )
#save(RMA_PNW_1989,file="/dmine/code/git/paper-PHD1/data/RMA_Rda/PNW_RMA_1989_2000.Rda")
write.csv(RMA_PNW_1989, file = "/dmine/code/git/paper-PHD1/data/RMA_csv/RMA_PNW_1989_2000.csv")
RMA_PNW$lossperacre <- RMA_PNW$loss / RMA_PNW$acres
RMA_PNW <- subset(RMA_PNW, month != "")
RMA_PNW$cropyear <- RMA_PNW$year
RMA_PNW$cropyear[RMA_PNW$monthcode >= 10] <- RMA_PNW$year + 1
RMA_PNW <- subset(RMA_PNW, year != 2016)
RMA_commodity_loss <- aggregate(RMA$loss, by = list(RMA$year, RMA$state, RMA$county, RMA$commodity), FUN = sum)
colnames(RMA_commodity_loss) <- c("year", "state", "county", "commodity", "loss")
RMA_commodity_loss <- subset(RMA_commodity_loss, commodity != "ADJUSTED GROSS REVENUE")
RMA_commodity_count <- aggregate(RMA$loss, by = list(RMA$year, RMA$state, RMA$county, RMA$commodity), FUN = length)
colnames(RMA_commodity_count) <- c("year", "state", "county", "commodity", "count")
RMA_commodity_count <- subset(RMA_commodity_count, commodity != "ADJUSTED GROSS REVENUE")
RMA_commodity_acres <- aggregate(RMA$acres, by = list(RMA$year, RMA$state, RMA$county, RMA$commodity), FUN = sum)
colnames(RMA_commodity_acres) <- c("year", "state", "county", "commodity", "acres")
RMA_commodity_acres <- subset(RMA_commodity_acres, commodity != "ADJUSTED GROSS REVENUE")
RMA_damage_loss <- aggregate(RMA$loss, by = list(RMA$year, RMA$state, RMA$county, RMA$commodity, RMA$damagecause), FUN = sum)
colnames(RMA_damage_loss) <- c("year", "state", "county", "commodity", "damagecause", "loss")
RMA_damage_count <- aggregate(RMA$loss, by = list(RMA$year, RMA$state, RMA$county, RMA$commodity, RMA$damagecause), FUN = length)
colnames(RMA_damage_count) <- c("year", "state", "county", "commodity", "damagecause", "count")
RMA_damage_acres <- aggregate(RMA$acres, by = list(RMA$year, RMA$state, RMA$county, RMA$commodity, RMA$damagecause), FUN = sum)
colnames(RMA_damage_acres) <- c("year", "state", "county", "commodity", "damaecause", "acres")
RMA_commodity_combined <- cbind(RMA_commodity_loss, RMA_commodity_count$count, RMA_commodity_acres$acres)
colnames(RMA_commodity_combined) <- c("year", "state", "county", "commodity", "loss", "count", "acres")
RMA_commodity_combined$lossperacre <- RMA_commodity_combined$loss / RMA_commodity_combined$acres
RMA_commodity_combined$lossperclaim <- RMA_commodity_combined$loss / RMA_commodity_combined$count
RMA_commodity_combined$acresperclaim <- RMA_commodity_combined$acres / RMA_commodity_combined$count
RMA_damage_combined <- cbind(RMA_damage_loss, RMA_damage_count$count, RMA_damage_acres$acres)
colnames(RMA_damage_combined) <- c("year", "state", "county", "commodity", "damagecause", "loss", "count", "acres")
RMA_damage_combined$lossperacre <- RMA_damage_combined$loss / RMA_damage_combined$acres
RMA_damage_combined$lossperclaim <- RMA_damage_combined$loss / RMA_damage_combined$count
RMA_damage_combined$acresperclaim <- RMA_damage_combined$acres / RMA_damage_combined$count
RMA_damage_combined$lossperacre[!is.finite(RMA_damage_combined$lossperacre)] <- 0
#write.csv(RMA_commodity_combined, file = paste("/nfs/soilsesfeedback-data/data/RMA/RMA_commodity_combined.csv", sep=""))
#write.csv(RMA_damage_combined, file = paste("/nfs/soilsesfeedback-data/data/RMA/RMA_damage_combined.csv", sep=""))
#subset for PNW - 2001 - 2015
RMA_damage_PNW <- subset(RMA_damage_combined, state == "WA" | state == "OR" | state == "ID")
RMA_damage_PNW <- subset(RMA_damage_PNW, county != "All Other Counties")
#write.csv(RMA_damage_PNW, file = "/dmine/data/USDA/crop_indemnity_PNW/RMA_damage_PNW_2001_2015.csv")
RMA_Idaho <- subset(RMA_damage_PNW, state == "ID")
RMA_Oregon <- subset(RMA_damage_PNW, state == "OR")
RMA_Washington <- subset(RMA_damage_PNW, state == "WA")
RMA_Idaho_IPNW <- subset(RMA_Idaho, county == "Idaho" | county == "Lewis" | county == "Nez Perce" | county == "Latah" | county == "Benewah")
RMA_Oregon_IPNW <- subset(RMA_Oregon, county == "Wasco" | county == "Sherman" | county == "Gilliam" | county == "Morrow" | county == "Umatilla" | county == "Union" | county == "Wallowa")
RMA_Washington_IPNW <- subset(RMA_Washington, county == "Douglas" | county == "Grant" | county == "Benton" | county == "Franklin" | county == "Walla Walla" | county == "Adams" | county == "Lincoln" | county == "Spokane" | county == "Whitman" | county == "Columbia" | county == "Garfield" | county == "Asotin")
RMA_damage_IPNW <- rbind(RMA_Idaho_IPNW,RMA_Oregon_IPNW,RMA_Washington_IPNW)
#-load summarized data
Southern_ID_sumloss <- read.csv(text=getURL
("https://raw.githubusercontent.com/erichseamon/dmine_anova/master/PNW_summary_all.csv"),
header = TRUE)
Southern_ID_sumloss_all_sum <- RMA_damage_PNW
Southern_ID_sumloss_all_sum <- aggregate(loss ~ year +
damagecause + state + county + commodity,
Southern_ID_sumloss, sum)
Southern_ID_count_all_count <- aggregate(count ~ year +
damagecause + county + commodity,
Southern_ID_sumloss, sum)
Southern_ID_sumloss_all_sum <-
Southern_ID_sumloss_all_sum[Southern_ID_sumloss_all_sum$loss >= 1, ]
#---subsetting for palouse/IPNW for sum, count, and lossperacre
#- SUM: sumloss_palouse_sum
#- COUNT: count_palouse_count
#- LOSSPERACRE: palouse_lossperacres
palousecounties <- c("Idaho", "Lewis", "Nez Perce", "Clearwater", "Latah", "Benewah", "Kootenai","Douglas", "Grant", "Benton", "Franklin", "Walla Walla", "Adams", "Lincoln", "Spokane", "Whitman", "Columbia", "Garfield", "Asotin","Wasco", "Sherman", "Gilliam", "Morrow", "Umatilla", "Union", "Wallowa")
sumloss_palouse <- Southern_ID_sumloss[Southern_ID_sumloss$county %in% palousecounties, ]
sumloss_nonpalouse <- Southern_ID_sumloss[!Southern_ID_sumloss$county %in% palousecounties, ]
sumloss_palouse_sum <- aggregate(loss ~ year +
damagecause + county + state + commodity,
sumloss_palouse, sum)
count_palouse_count <- aggregate(count ~ year +
damagecause + county + state + commodity,
sumloss_palouse, sum)
palouse_lossperacres <- aggregate(lossperacres ~ year +
damagecause + county + state + commodity,
sumloss_palouse, mean)
sumloss_palouse_sum <-
sumloss_palouse_sum[sumloss_palouse_sum$loss >= 1, ]
#creating transforms for sum loss - cube root, square root, inverse, log10 and natural log, and then scale and center
Math.cbrt <- function(x) {
sign(x) * abs(x)^(1/3)
}
#-calculate cube loss
sumloss_palouse_sum$cube_loss <- Math.cbrt(sumloss_palouse_sum$loss)
#-remove zeros
sumloss_palouse_sum <- subset(sumloss_palouse_sum, loss > 0)
#-use a log transform
sumloss_palouse_sum$log10_loss <- log10(sumloss_palouse_sum$loss)
sumloss_palouse_sum$log_loss <- log(sumloss_palouse_sum$loss)
#-inverse transform
sumloss_palouse_sum$inverse_loss <- 1/sumloss_palouse_sum$loss
#sq root transformation
sumloss_palouse_sum$sqroot_loss <- sqrt(sumloss_palouse_sum$loss)
#--scale and center
sumloss_palouse_sum$scaled_cube_loss <- scale(sumloss_palouse_sum$cube_loss)
sumloss_palouse_sum$scaled_inverse_loss <-
scale(sumloss_palouse_sum$inverse_loss, center = TRUE, scale = TRUE)
#----
#-Loading all WHEAT claims for the palouse from 1989-2015
palouse_sumloss <- read.csv(text=getURL
("https://raw.githubusercontent.com/erichseamon/dmine_anova/master/Palouse_summary_sumloss.csv"),
header = TRUE)
palouse_counts <- read.csv(text=getURL
("https://raw.githubusercontent.com/erichseamon/dmine_anova/master/Palouse_summary_counts.csv"),
header = TRUE)
#use a cube transformation on loss for WHEAT claims
Math.cbrt <- function(x) {
sign(x) * abs(x)^(1/3)
}
Southern_ID_sumloss_all_sum$cube_loss <- Math.cbrt(Southern_ID_sumloss_all_sum$loss)
Southern_ID_count_all_count$cube_counts <- Math.cbrt(Southern_ID_count_all_count$count)
#--aggregate palouse but not by county
palouse_sumloss_aggregate_county <- aggregate(loss ~ damagecause + year + commodity + county,
palouse_sumloss, mean)
#--aggregate palouse
palouse_sumloss_aggregate <- aggregate(loss ~ damagecause + year + commodity + county,
palouse_sumloss, sum)
#-calculate cube loss
palouse_sumloss_aggregate$cube_loss <- Math.cbrt(palouse_sumloss_aggregate$loss)
#-remove zeros
palouse_sumloss_aggregate <- subset(palouse_sumloss_aggregate, loss > 0)
#-use a log transform
palouse_sumloss_aggregate$log10_loss <- log10(palouse_sumloss_aggregate$loss)
palouse_sumloss_aggregate$log_loss <- log(palouse_sumloss_aggregate$loss)
#-inverse transform
palouse_sumloss_aggregate$inverse_loss <- 1/palouse_sumloss_aggregate$loss
#sq root transformation
palouse_sumloss_aggregate$sqroot_loss <- sqrt(palouse_sumloss_aggregate$loss)
#--scale and center
palouse_sumloss_aggregate$scaled_cube_loss <- scale(palouse_sumloss_aggregate$cube_loss)
palouse_sumloss_aggregate$scaled_inverse_loss <-
scale(palouse_sumloss_aggregate$inverse_loss, center = TRUE, scale = TRUE)
#----Step 2. Original Insurance Loss Dataset - Pacific Northwest
In Step 2, we document the following table, which is an example of the original data that was acquired from the USDA Risk Management Agency (RMA). Each record represents an individual insurance claim.
knitr::kable(RMA_PNW[1:10,], format="markdown") %>%
kable_styling(bootstrap_options = c("striped", "hover", "condensed"))| year | state | county | commodity | damagecause | monthcode | month | acres | loss | lossperacre | cropyear | |
|---|---|---|---|---|---|---|---|---|---|---|---|
| 11338 | 2001 | ID | Ada | All Other Crops | Drought | 9 | SEP | 17.000 | 153.00 | 9.000000 | 2001 |
| 11339 | 2001 | ID | Ada | All Other Crops | Heat | 8 | AUG | 105.200 | 5249.00 | 49.895437 | 2001 |
| 11340 | 2001 | ID | Ada | All Other Crops | Freeze | 4 | APR | 125.000 | 4500.00 | 36.000000 | 2001 |
| 11341 | 2001 | ID | Ada | All Other Crops | Wind/Excess Wind | 5 | MAY | 50.000 | 1800.00 | 36.000000 | 2001 |
| 11342 | 2001 | ID | Ada | All Other Crops | Wind/Excess Wind | 4 | APR | 92.500 | 3330.00 | 36.000000 | 2001 |
| 11343 | 2001 | ID | Bannock | WHEAT | Drought | 8 | AUG | 133.000 | 1212.00 | 9.112782 | 2001 |
| 11344 | 2001 | ID | Bannock | WHEAT | Drought | 9 | SEP | 777.520 | 24807.00 | 31.905289 | 2001 |
| 11345 | 2001 | ID | Bannock | WHEAT | Drought | 7 | JUL | 3529.754 | 54726.46 | 15.504327 | 2001 |
| 11346 | 2001 | ID | Bannock | WHEAT | Heat | 7 | JUL | 19.796 | 2371.60 | 119.801980 | 2001 |
| 11347 | 2001 | ID | Bannock | WHEAT | Heat | 8 | AUG | 25.000 | 904.00 | 36.160000 | 2001 |
Step 3. Aggregated Insurance Loss Dataset - Pacific Northwest
In Step 3, the following table is an example of an aggregated dataset derived from the original insurance loss csv files. Here we have summarized claims by year, county, commodity, and damage cause. Each unique combination is summarized, which includes the total summarized loss, the number of claims, the total summarized acreage, loss per acre, loss per claim, and acres per claim. This dataset was the basis for our data examination.
knitr::kable(RMA_damage_PNW[1:10,], format="markdown") %>%
kable_styling(bootstrap_options = c("striped", "hover", "condensed"))| year | state | county | commodity | damagecause | loss | count | acres | lossperacre | lossperclaim | acresperclaim | |
|---|---|---|---|---|---|---|---|---|---|---|---|
| 11 | 2014 | ID | Ada | All Other Crops | Area Plan Crops Only | 1398 | 1 | 0 | 0.0000000 | 1398.0 | 0 |
| 12 | 2015 | ID | Ada | All Other Crops | Area Plan Crops Only | 11810 | 1 | 0 | 0.0000000 | 11810.0 | 0 |
| 211 | 2008 | OR | Baker | All Other Crops | Area Plan Crops Only | 15292 | 2 | 5482 | 2.7894929 | 7646.0 | 2741 |
| 212 | 2010 | OR | Baker | All Other Crops | Area Plan Crops Only | 1819 | 2 | 2282 | 0.7971078 | 909.5 | 1141 |
| 215 | 2009 | ID | Bannock | All Other Crops | Area Plan Crops Only | 2284 | 1 | 600 | 3.8066667 | 2284.0 | 600 |
| 245 | 2008 | ID | Bear Lake | All Other Crops | Area Plan Crops Only | 8102 | 1 | 2468 | 3.2828201 | 8102.0 | 2468 |
| 246 | 2012 | ID | Bear Lake | All Other Crops | Area Plan Crops Only | 7459 | 1 | 2468 | 3.0222853 | 7459.0 | 2468 |
| 291 | 2010 | ID | Bingham | All Other Crops | Area Plan Crops Only | 4197 | 1 | 192 | 21.8593750 | 4197.0 | 192 |
| 292 | 2011 | ID | Bingham | All Other Crops | Area Plan Crops Only | 7769 | 1 | 240 | 32.3708333 | 7769.0 | 240 |
| 293 | 2012 | ID | Bingham | All Other Crops | Area Plan Crops Only | 6786 | 1 | 168 | 40.3928571 | 6786.0 | 168 |
Step 4: inland Pacific Northwest (iPNW) Study Area
In Step 4, we document the 24 county iPNW study area.
library(data.table)
library(maptools)
library(classInt)
library(leaflet)
library(dplyr)
library(raster)
library(htmltools)
states <- readShapePoly('/dmine/data/states/states.shp',
proj4string=CRS
("+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0"))
projection = CRS("+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0")
#setwd("/dmine/data/counties/")
counties <- readShapePoly('/dmine/data/counties/threestate_palouse.shp',
proj4string=CRS
("+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0"))
projection = CRS("+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0")
#pal <- colorNumeric(palette = c("white", "orange", "darkorange", "red", "darkred"),
# domain = counties$NAME)
exte <- as.vector(extent(counties))
#label <- paste(sep = "<br/>", counties$NAME, round(counties$NAME, 0))
#markers <- data.frame(label)
#labs <- as.list(counties$NAME)
counties <- counties[counties$NAME != "Kootenai",]
counties <- counties[counties$NAME != "Clearwater",]
lat_long <- coordinates(counties)
counties_a <- data.frame(counties)
labs <- lapply(seq(nrow(counties_a)), function(i) {
paste0(as.character(counties_a[i, "NAME"]))
})
#leaflet(data = counties, options = leafletOptions(minZoom = 6, maxZoom = 6)) %>% addProviderTiles("Stamen.TonerLite") %>% fitBounds(exte[1], exte[3], exte[2], exte[4]) %>% addPolygons(color = "blue", weight = 1) %>%
#--
mapz <- leaflet(data = counties, options = leafletOptions(zoomControl = FALSE, zoomSnap = .4, zoomDelta = .4)) %>%
addProviderTiles(providers$Hydda.Base) %>%
addProviderTiles(providers$Stamen.TonerLines) %>%
fitBounds(exte[1], exte[3], exte[2], exte[4]) %>% addPolygons(color = "lightyellow", fillOpacity = .8, weight = 1) %>% addPolygons(data = states, fillOpacity = 0, weight = 4, stroke = TRUE, smoothFactor = 0.5, color = "black") %>% addPolygons(color = "gray", fillOpacity = 0, weight = 2) %>%
# addMiniMap(
# tiles = providers$Stamen.TonerLite,
# position = 'topleft',
# width = 175, height = 175,
# toggleDisplay = FALSE) %>%
addLabelOnlyMarkers(data = counties, lng = lat_long[,1], lat = lat_long[,2], label = lapply(labs, HTML), labelOptions = labelOptions(style = list("font-family" = "serif"), noHide = TRUE, direction = 'middle', textOnly = TRUE, textsize = "16px", color = "white", offset = c(20, 0), markerOptions(riseOnHover = TRUE, col = "white")
))
addScaleBar(mapz, position = c("topright"), options = scaleBarOptions())#markerOptions(map, riseOnHover = TRUE)
# addLegend(colors = ~NAME, labels = ~NAME, values = ~NAME, opacity = 0.5, title = NULL,
# position = "bottomright")Step 5: PCA: PNW insurance loss, by COUNTY with COMMODITY loadings: 2001 to 2015
In Step 5, we perform a PCA for the insurance loss for the entire PNW, by county, with commodity as the factor loadings. Data from 2001 is 2015 is used. We additionally have generated a scree plot that shows the proportion of variance explained by the individual components.
#PNW by commodity
library(reshape2) ; RMA_damage_PNW_transformed_commodity <- dcast(RMA_damage_PNW, county + state ~ commodity, value.var = c("loss"), sum)
RMA_damage_PNW_transformed_commodity$statecounty <- paste(RMA_damage_PNW_transformed_commodity$county, "_", RMA_damage_PNW_transformed_commodity$state, sep="")
rownames(RMA_damage_PNW_transformed_commodity) <- RMA_damage_PNW_transformed_commodity$statecounty
RMA_damage_PNW_exogenous_commodity <- RMA_damage_PNW_transformed_commodity[,3:41]
RMA_damage_PNW_exogenous_commodity <- RMA_damage_PNW_exogenous_commodity[-1]
RMA_damage_PNW_exogenous_commodity <- RMA_damage_PNW_exogenous_commodity[-20]
RMA_damage_PNW_exogenous_commodity <- RMA_damage_PNW_exogenous_commodity[-37]
fit <- princomp(RMA_damage_PNW_exogenous_commodity, cor=TRUE)
#autoplot(prcomp(RMA_damage_PNW_exogenous_commodity, scale = TRUE, center = TRUE), data = RMA_damage_PNW_transformed_commodity, colour = 'county', loadings = TRUE, loadings.label = TRUE, loadings.label.size = 3, label = FALSE, label.size = 3, legend = FALSE) + geom_text(vjust=-1, label=rownames(RMA_damage_PNW_exogenous_commodity)) + theme(legend.position = "none")
autoplot(prcomp(RMA_damage_PNW_exogenous_commodity, center = TRUE, scale = TRUE), frame = TRUE, frame.type = 'norm', data = RMA_damage_PNW_transformed_commodity, colour = 'state', loadings = TRUE, loadings.label = TRUE, loadings.label.size = 3, label = FALSE, label.size = 2, legend = FALSE) + theme_bw() + geom_text(vjust=-1, label=rownames(RMA_damage_PNW_exogenous_commodity)) + theme(legend.position = "none") + ggtitle("Principle Components: PNW insurance loss\nby county with commodity loadings ") + theme(plot.title = element_text(size =28, family = "Serif")) + theme(axis.title.y = element_text(family = "Serif", size=16), axis.title.x = element_text(family = "Serif", size = 16), axis.text.x = element_text(size=rel(1.5), hjust = 1, family = "Serif"), axis.text.y = element_text(size=rel(1.5), hjust = 1, family = "Serif"))
#pca_RMA_damage_IPNW_exogenous$rotation
#plot(pca_RMA_damage_IPNW_exogenous$rotation)
std_dev <- fit$sdev
pr_var <- std_dev^2
#pr_var[1:10]
prop_varex <- pr_var/sum(pr_var)
plot(prop_varex, xlab = "Principal Component",
ylab = "Proportion of Variance Explained",
type = "b", main = "Scree Plot")
Step 6: PCA: PNW insurance loss, by YEAR with DAMAGE CAUSE loadings: 2001 to 2015
In Step 6, we perform a PCA for the insurance loss for the entire PNW, for all commodities by year, with damage cause as the factor loadings. Data from 2001 is 2015 is used. We additionally have generated a scree plot that shows the proportion of variance explained by the individual components.
library(reshape2) ; RMA_damage_PNW_transformed_year <- dcast(RMA_damage_PNW, year ~ damagecause, value.var = c("loss"), sum)
#RMA_damage_PNW_transformed_year$statecounty <- paste(RMA_damage_PNW_transformed_year$county, "_", RMA_damage_PNW_transformed_year$state, sep="")
rownames(RMA_damage_PNW_transformed_year) <- RMA_damage_PNW_transformed_year$year
RMA_damage_PNW_exogenous_year <- RMA_damage_PNW_transformed_year[,2:31]
#PNW using year as group - with damage cause as loadings
#fit <- princomp(RMA_damage_PNW_exogenous_year, cor=TRUE, scale = TRUE, center = TRUE)
RMA_damage_PNW_exogenous_year_factor <- cbind(RMA_damage_PNW_exogenous_year[,2:3], RMA_damage_PNW_exogenous_year[6], RMA_damage_PNW_exogenous_year[8], RMA_damage_PNW_exogenous_year[,11:14], RMA_damage_PNW_exogenous_year[16:17])
fit <- princomp(RMA_damage_PNW_exogenous_year_factor, cor=TRUE)
autoplot(prcomp(RMA_damage_PNW_exogenous_year_factor, center = TRUE, scale = TRUE), data = RMA_damage_PNW_exogenous_year_factor, colour = rownames(RMA_damage_PNW_exogenous_year_factor), loadings = TRUE, loadings.label = TRUE, loadings.label.size = 3, label = FALSE, label.size = 2, legend = FALSE) + theme_bw() + geom_text(vjust=-1, label=rownames(RMA_damage_PNW_exogenous_year_factor)) + theme(legend.position = "none") + ggtitle("Principle Components: PNW insurance loss\nby year with damage cause loadings ") + theme(plot.title = element_text(size =28, family = "Serif")) + theme(axis.title.y = element_text(family = "Serif", size=16), axis.title.x = element_text(family = "Serif", size = 16), axis.text.x = element_text(size=rel(1.5), hjust = 1, family = "Serif"), axis.text.y = element_text(size=rel(1.5), hjust = 1, family = "Serif"))
#pca_RMA_damage_IPNW_exogenous$rotation
#plot(pca_RMA_damage_IPNW_exogenous$rotation)
std_dev <- fit$sdev
pr_var <- std_dev^2
#pr_var[1:10]
prop_varex <- pr_var/sum(pr_var)
plot(prop_varex, xlab = "Principal Component",
ylab = "Proportion of Variance Explained",
type = "b", main = "Scree Plot")
Step 7: PCA: PNW insurance loss, by MONTH with DAMAGE CAUSE loadings: 2001 to 2015
In Step 7, we perform a PCA for the insurance loss for the entire PNW, for all commodities by month, with damage cause as the factor loadings. Data from 2001 is 2015 is used. We additionally have generated a scree plot that shows the proportion of variance explained by the individual components.
#PNW damage cause as loadings using months as groups
#remove negatives related to lossperacre calc
RMA_PNW_noneg <- subset(RMA_PNW, lossperacre != Inf)
RMA_PNW_noneg2 <- subset(RMA_PNW_noneg, lossperacre >=0 )
RMA_PNW_noneg3 <- subset(RMA_PNW_noneg2, month !="" )
library(reshape2) ; RMA_damage_PNW_transformed_month <- dcast(RMA_PNW_noneg3, month ~ damagecause, value.var = c("loss"), mean)
is.nan.data.frame <- function(x)
do.call(cbind, lapply(x, is.nan))
RMA_damage_PNW_transformed_month[is.nan(RMA_damage_PNW_transformed_month)] <- 0
#RMA_damage_PNW_transformed_year$statecounty <- paste(RMA_damage_PNW_transformed_year$county, "_", RMA_damage_PNW_transformed_year$state, sep="")
rownames(RMA_damage_PNW_transformed_month) <- RMA_damage_PNW_transformed_month$month
RMA_damage_PNW_exogenous_month <- RMA_damage_PNW_transformed_month[,2:31]
RMA_damage_PNW_exogenous_month <- cbind(RMA_damage_PNW_exogenous_month[,2:3], RMA_damage_PNW_exogenous_month[6], RMA_damage_PNW_exogenous_month[8], RMA_damage_PNW_exogenous_month[,11:14], RMA_damage_PNW_exogenous_month[16:17])
#autoplot(prcomp(RMA_damage_PNW_exogenous_year, center = TRUE, scale = TRUE), data = RMA_damage_PNW_transformed_year, colour = 'year', loadings = TRUE, loadings.label = TRUE, loadings.label.size = 3, label = FALSE, label.size = 3, legend = FALSE) + geom_text(vjust=-1, label=rownames(RMA_damage_PNW_exogenous_year)) + theme(legend.position = "none")
#PNW using month as group - with damage cause as loadings
fit <- princomp(RMA_damage_PNW_exogenous_month, cor=TRUE)
autoplot(prcomp(RMA_damage_PNW_exogenous_month, center = TRUE, scale = TRUE), data = RMA_damage_PNW_transformed_month, colour = 'month', loadings = TRUE, loadings.label = TRUE, loadings.label.size = 3, label = FALSE, label.size = 2, legend = FALSE) + theme_bw() + geom_text(vjust=-1, label=rownames(RMA_damage_PNW_exogenous_month)) + theme(legend.position = "none") + ggtitle("Principle Components: PNW insurance loss\nby month with damage cause loadings ") + theme(plot.title = element_text(size =28, family = "Serif")) + theme(axis.title.y = element_text(family = "Serif", size=16), axis.title.x = element_text(family = "Serif", size = 16), axis.text.x = element_text(size=rel(1.5), hjust = 1, family = "Serif"), axis.text.y = element_text(size=rel(1.5), hjust = 1, family = "Serif"))
#pca_RMA_damage_IPNW_exogenous$rotation
#plot(pca_RMA_damage_IPNW_exogenous$rotation)
std_dev <- fit$sdev
pr_var <- std_dev^2
#pr_var[1:10]
prop_varex <- pr_var/sum(pr_var)
plot(prop_varex, xlab = "Principal Component",
ylab = "Proportion of Variance Explained",
type = "b", main = "Scree Plot")
Step 8: PCA: PNW WHEAT insurance loss, by YEAR with DAMAGE CAUSE loadings: 2001 to 2015
In Step 8, we perform a PCA for the insurance loss for the entire PNW, for wheat by year, with damage cause as the factor loadings. Data from 2001 is 2015 is used. We additionally have generated a scree plot that shows the proportion of variance explained by the individual components.
#PNW damage cause as loadings using years as groups, JUST FOR WHEAT
RMA_damage_PNW_wheat_year <- subset(RMA_damage_PNW, commodity == "WHEAT")
library(reshape2) ; RMA_damage_PNW_transformed_wheat_year <- dcast(RMA_damage_PNW_wheat_year, year ~ damagecause, value.var = c("loss"), sum)
#RMA_damage_PNW_transformed_year$statecounty <- paste(RMA_damage_PNW_transformed_year$county, "_", RMA_damage_PNW_transformed_year$state, sep="")
rownames(RMA_damage_PNW_transformed_wheat_year) <- RMA_damage_PNW_transformed_wheat_year$year
RMA_damage_PNW_exogenous_wheat_year <- RMA_damage_PNW_transformed_wheat_year
par(mar=c(8, 6.1, 4, 2.1), family = 'serif', mgp=c(4, 1, 0), las=0)
#PNW using year as group - with damage cause as loadings: JUST FOR WHEAT
RMA_damage_PNW_exogenous_wheat_year_factor <- cbind(RMA_damage_PNW_exogenous_wheat_year[,3:4], RMA_damage_PNW_exogenous_wheat_year[7], RMA_damage_PNW_exogenous_wheat_year[9], RMA_damage_PNW_exogenous_wheat_year[,12:15], RMA_damage_PNW_exogenous_wheat_year[17:18])
fit <- princomp(RMA_damage_PNW_exogenous_wheat_year_factor, cor=TRUE)
autoplot(prcomp(RMA_damage_PNW_exogenous_wheat_year_factor, center = TRUE, scale = TRUE), data = RMA_damage_PNW_exogenous_wheat_year_factor, colour = rownames(RMA_damage_PNW_exogenous_wheat_year_factor), loadings = TRUE, loadings.label = TRUE, loadings.label.size = 3, label = FALSE, label.size = 2, legend = FALSE) + theme_bw() + geom_text(vjust=-1, label=rownames(RMA_damage_PNW_exogenous_wheat_year_factor)) + theme(legend.position = "none") + ggtitle("Principle Components: PNW WHEAT insurance loss\nby year with damage cause loadings ") + theme(plot.title = element_text(size =28, family = "Serif")) + theme(axis.title.y = element_text(family = "Serif", size=16), axis.title.x = element_text(family = "Serif", size = 16), axis.text.x = element_text(size=rel(1.5), hjust = 1, family = "Serif"), axis.text.y = element_text(size=rel(1.5), hjust = 1, family = "Serif"))
#pca_RMA_damage_IPNW_exogenous$rotation
#plot(pca_RMA_damage_IPNW_exogenous$rotation)
std_dev <- fit$sdev
pr_var <- std_dev^2
#pr_var[1:10]
prop_varex <- pr_var/sum(pr_var)
plot(prop_varex, xlab = "Principal Component",
ylab = "Proportion of Variance Explained",
type = "b", main = "Scree Plot")
Step 9: PCA: PNW WHEAT insurance loss, by COUNTY with DAMAGE CAUSE loadings: 2001 to 2015
In Step 9, we perform a PCA for the insurance loss for the entire PNW, for wheat by county, with damage cause as the factor loadings. Data from 2001 is 2015 is used. We additionally have generated a scree plot that shows the proportion of variance explained by the individual components.
#PNW wheat
RMA_damage_PNW_wheat <- subset(RMA_damage_PNW, commodity == "WHEAT")
library(reshape2) ; RMA_damage_PNW_transformed_wheat <- dcast(RMA_damage_PNW_wheat, county + state ~ damagecause, value.var = c("loss"), sum)
RMA_damage_PNW_transformed_wheat$statecounty <- paste(RMA_damage_PNW_transformed_wheat$county, "_", RMA_damage_PNW_transformed_wheat$state, sep="")
rownames(RMA_damage_PNW_transformed_wheat) <- RMA_damage_PNW_transformed_wheat$statecounty
#RMA_damage_PNW_exogenous_wheat <- RMA_damage_PNW_transformed_wheat[,3:27]
RMA_damage_PNW_exogenous_wheat <- RMA_damage_PNW_transformed_wheat
#autoplot(prcomp(RMA_damage_PNW_exogenous_month), data = RMA_damage_PNW_transformed_month, colour = 'month', loadings = TRUE, loadings.label = TRUE, loadings.label.size = 3, label = FALSE, label.size = 3, legend = FALSE) + geom_text(vjust=-1, label=rownames(RMA_damage_PNW_exogenous_month)) + theme(legend.position = "none")
#PNW wheat
RMA_damage_PNW_exogenous_wheat <- cbind(RMA_damage_PNW_exogenous_wheat[,4:5], RMA_damage_PNW_exogenous_wheat[8], RMA_damage_PNW_exogenous_wheat[10], RMA_damage_PNW_exogenous_wheat[,13:16], RMA_damage_PNW_exogenous_wheat[18:19])
fit <- princomp(RMA_damage_PNW_exogenous_wheat, cor=TRUE)
#autoplot(prcomp(RMA_damage_PNW_exogenous_wheat, scale = TRUE, center = TRUE), data = RMA_damage_PNW_transformed_wheat, colour = 'county', loadings = TRUE, loadings.label = TRUE, loadings.label.size = 3, label = FALSE, label.size = 3, legend = FALSE) + geom_text(vjust=-1, label=rownames(RMA_damage_PNW_exogenous_wheat)) + theme(legend.position = "none")
autoplot(prcomp(RMA_damage_PNW_exogenous_wheat, center = TRUE, scale = TRUE), frame = TRUE, frame.type = 'norm', data = RMA_damage_PNW_transformed_wheat, colour = 'state', loadings = TRUE, loadings.label = TRUE, loadings.label.size = 3, label = FALSE, label.size = 2, legend = FALSE) + theme_bw() + geom_text(vjust=-1, label=rownames(RMA_damage_PNW_exogenous_wheat)) + theme(legend.position = "none") + ggtitle("Principle Components: PNW wheat insurance loss\nby county with damage cause loadings ") + theme(plot.title = element_text(size =28, family = "Serif")) + theme(axis.title.y = element_text(family = "Serif", size=16), axis.title.x = element_text(family = "Serif", size = 16), axis.text.x = element_text(size=rel(1.5), hjust = 1, family = "Serif"), axis.text.y = element_text(size=rel(1.5), hjust = 1, family = "Serif"))
#pca_RMA_damage_IPNW_exogenous$rotation
#plot(pca_RMA_damage_IPNW_exogenous$rotation)
std_dev <- fit$sdev
pr_var <- std_dev^2
#pr_var[1:10]
prop_varex <- pr_var/sum(pr_var)
plot(prop_varex, xlab = "Principal Component",
ylab = "Proportion of Variance Explained",
type = "b", main = "Scree Plot")
Step 10: PCA: PNW APPLES insurance loss, by COUNTY with DAMAGE CAUSE loadings: 2001 to 2015
In Step 5, we perform a PCA for the insurance loss for the entire PNW, for apples by county, with damage cause as the factor loadings. Data from 2001 is 2015 is used. We additionally have generated a scree plot that shows the proportion of variance explained by the individual components.
#PNW apples
RMA_damage_PNW_apples <- subset(RMA_damage_PNW, commodity == "APPLES")
library(reshape2) ; RMA_damage_PNW_transformed_apples <- dcast(RMA_damage_PNW_apples, county + state ~ damagecause, value.var = c("loss"), sum)
RMA_damage_PNW_transformed_apples$statecounty <- paste(RMA_damage_PNW_transformed_apples$county, "_", RMA_damage_PNW_transformed_apples$state, sep="")
rownames(RMA_damage_PNW_transformed_apples) <- RMA_damage_PNW_transformed_apples$statecounty
RMA_damage_PNW_exogenous_apples <- RMA_damage_PNW_transformed_apples[,3:24]
RMA_damage_PNW_exogenous_apples <- cbind(RMA_damage_PNW_exogenous_apples[,1:4], RMA_damage_PNW_exogenous_apples[,7:10], RMA_damage_PNW_exogenous_apples[,13:14])
fit <- princomp(RMA_damage_PNW_exogenous_apples, cor=TRUE)
#PNW apples
#autoplot(prcomp(RMA_damage_PNW_exogenous_apples), data = RMA_damage_PNW_transformed_apples, colour = 'county', loadings = TRUE, loadings.label = TRUE, loadings.label.size = 3, label = FALSE, label.size = 3, legend = FALSE) + geom_text(vjust=-1, label=rownames(RMA_damage_PNW_exogenous_apples)) + theme(legend.position = "none")
autoplot(prcomp(RMA_damage_PNW_exogenous_apples, center = TRUE, scale = TRUE), frame = TRUE, frame.type = 'norm', data = RMA_damage_PNW_transformed_apples, colour = 'state', loadings = TRUE, loadings.label = TRUE, loadings.label.size = 3, label = FALSE, label.size = 2, legend = FALSE) + theme_bw() + geom_text(vjust=-1, label=rownames(RMA_damage_PNW_exogenous_apples)) + theme(legend.position = "none") + ggtitle("Principle Components: PNW apples insurance loss\nby county with damage cause loadings ") + theme(plot.title = element_text(size =28, family = "Serif")) + theme(axis.title.y = element_text(family = "Serif", size=16), axis.title.x = element_text(family = "Serif", size = 16), axis.text.x = element_text(size=rel(1.5), hjust = 1, family = "Serif"), axis.text.y = element_text(size=rel(1.5), hjust = 1, family = "Serif"))
#pca_RMA_damage_IPNW_exogenous$rotation
#plot(pca_RMA_damage_IPNW_exogenous$rotation)
std_dev <- fit$sdev
pr_var <- std_dev^2
#pr_var[1:10]
prop_varex <- pr_var/sum(pr_var)
plot(prop_varex, xlab = "Principal Component",
ylab = "Proportion of Variance Explained",
type = "b", main = "Scree Plot")
Step 11: PCA: PNW BARLEY insurance loss, by COUNTY with DAMAGE CAUSE loadings: 2001 to 2015
In Step 11, we perform a PCA for the insurance loss for the entire PNW, for barley by county, with damage cause as the factor loadings. Data from 2001 is 2015 is used. We additionally have generated a scree plot that shows the proportion of variance explained by the individual components.
#PNW barley
RMA_damage_PNW_barley <- subset(RMA_damage_PNW, commodity == "BARLEY")
library(reshape2) ; RMA_damage_PNW_transformed_barley <- dcast(RMA_damage_PNW_barley, county + state ~ damagecause, value.var = c("loss"), sum)
RMA_damage_PNW_transformed_barley$statecounty <- paste(RMA_damage_PNW_transformed_barley$county, "_", RMA_damage_PNW_transformed_barley$state, sep="")
rownames(RMA_damage_PNW_transformed_barley) <- RMA_damage_PNW_transformed_barley$statecounty
#RMA_damage_PNW_exogenous_barley <- RMA_damage_PNW_transformed_barley[,3:23]
RMA_damage_PNW_exogenous_barley <- RMA_damage_PNW_transformed_barley
RMA_damage_PNW_exogenous_barley <- cbind(RMA_damage_PNW_exogenous_barley[,3:4], RMA_damage_PNW_exogenous_barley[7], RMA_damage_PNW_exogenous_barley[9], RMA_damage_PNW_exogenous_barley[,12:15], RMA_damage_PNW_exogenous_barley[17:18])
fit <- princomp(RMA_damage_PNW_exogenous_barley, cor=TRUE)
#PNW barley
#autoplot(prcomp(RMA_damage_PNW_exogenous_barley), data = RMA_damage_PNW_transformed_barley, colour = 'county', loadings = TRUE, loadings.label = TRUE, loadings.label.size = 3, label = FALSE, label.size = 3, legend = FALSE) + geom_text(vjust=-1, label=rownames(RMA_damage_PNW_exogenous_barley)) + theme(legend.position = "none")
autoplot(prcomp(RMA_damage_PNW_exogenous_barley, center = TRUE, scale = TRUE), frame = TRUE, frame.type = 'norm', data = RMA_damage_PNW_transformed_barley, colour = 'state', loadings = TRUE, loadings.label = TRUE, loadings.label.size = 3, label = FALSE, label.size = 2, legend = FALSE) + theme_bw() + geom_text(vjust=-1, label=rownames(RMA_damage_PNW_exogenous_barley)) + theme(legend.position = "none") + ggtitle("Principle Components: PNW barley insurance loss\nby county with damage cause loadings ") + theme(plot.title = element_text(size =28, family = "Serif")) + theme(axis.title.y = element_text(family = "Serif", size=16), axis.title.x = element_text(family = "Serif", size = 16), axis.text.x = element_text(size=rel(1.5), hjust = 1, family = "Serif"), axis.text.y = element_text(size=rel(1.5), hjust = 1, family = "Serif"))
#pca_RMA_damage_IPNW_exogenous$rotation
#plot(pca_RMA_damage_IPNW_exogenous$rotation)
std_dev <- fit$sdev
pr_var <- std_dev^2
#pr_var[1:10]
prop_varex <- pr_var/sum(pr_var)
plot(prop_varex, xlab = "Principal Component",
ylab = "Proportion of Variance Explained",
type = "b", main = "Scree Plot")
Step 12: PCA: IPNW insurance loss by COUNTY with DAMAGE CAUSE loadings: 2001 to 2015
In Step 5, we perform a PCA for the insurance loss for the inland PNW, for all commodities by county, with damage cause as the factor loadings. Data from 2001 is 2015 is used. We additionally have generated a scree plot that shows the proportion of variance explained by the individual components.
#all PNW by damage
library(reshape2) ; RMA_damage_PNW_transformed <- dcast(RMA_damage_PNW, county + state ~ damagecause, value.var = c("loss"), sum)
RMA_damage_PNW_transformed$statecounty <- paste(RMA_damage_PNW_transformed$county, "_", RMA_damage_PNW_transformed$state, sep="")
rownames(RMA_damage_PNW_transformed) <- RMA_damage_PNW_transformed$statecounty
RMA_damage_PNW_exogenous <- RMA_damage_PNW_transformed[,3:32]
par(mar=c(8, 6.1, 4, 2.1), family = 'serif', mgp=c(4, 1, 0), las=0)
fit <- stats::princomp(RMA_damage_PNW_exogenous, cor=TRUE)
autoplot(prcomp(RMA_damage_PNW_exogenous, center = TRUE, scale = TRUE), frame = TRUE, frame.type = 'norm', data = RMA_damage_PNW_transformed, colour = 'state', loadings = TRUE, loadings.label = TRUE, loadings.label.size = 3, label = FALSE, label.size = 2, legend = FALSE) + theme_bw() + geom_text(vjust=-1, label=rownames(RMA_damage_PNW_exogenous)) + theme(legend.position = "none") + ggtitle("Principle Components: IPNW insurance loss\nby county with damage cause loadings: 2001 to 2015") + theme(plot.title = element_text(size =28, family = "Serif")) + theme(axis.title.y = element_text(family = "Serif", size=16), axis.title.x = element_text(family = "Serif", size = 16), axis.text.x = element_text(size=rel(1.5), hjust = 1, family = "Serif"), axis.text.y = element_text(size=rel(1.5), hjust = 1, family = "Serif"))
par(mar=c(8, 6.1, 4, 2.1), family = 'serif', mgp=c(4, 1, 0), las=0)#pca_RMA_damage_IPNW_exogenous$rotation
#plot(pca_RMA_damage_IPNW_exogenous$rotation)
std_dev <- fit$sdev
pr_var <- std_dev^2
#pr_var[1:10]
prop_varex <- pr_var/sum(pr_var)
plot(prop_varex, xlab = "Principal Component",
ylab = "Proportion of Variance Explained",
type = "b", main = "Scree Plot")
Step 13: PCA: IPNW insurance loss, by COUNTY with DAMAGE CAUSE loadings: 2001 to 2015
In Step 5, we perform a PCA for the insurance loss for the inland PNW, for all commodities by county, with damage cause as the factor loadings. Data from 2001 is 2015 is used. We additionally have generated a scree plot that shows the proportion of variance explained by the individual components.
#just the IPNW
RMA_damage_IPNW$log10loss <- log10(RMA_damage_IPNW$loss)
RMA_damage_IPNW_wheat <- subset(RMA_damage_IPNW, commodity == "WHEAT")
library(reshape2) ; RMA_damage_IPNW_transformed <- dcast(RMA_damage_IPNW_wheat, county + state ~ damagecause, value.var = c("loss"), sum)
RMA_damage_IPNW_transformed$statecounty <- paste(RMA_damage_IPNW_transformed$county, "_", RMA_damage_IPNW_transformed$state, sep="")
rownames(RMA_damage_IPNW_transformed) <- RMA_damage_IPNW_transformed$statecounty
RMA_damage_IPNW_exogenous <- RMA_damage_IPNW_transformed[,3:26]
#IPNW
#RMA_damage_IPNW_exogenous <- scale(RMA_damage_IPNW_exogenous, scale = TRUE, center = TRUE)
pca_RMA_damage_IPNW_exogenous <- prcomp(RMA_damage_IPNW_exogenous, scale = TRUE, center = TRUE)
RMA_damage_IPNW_exogenous_factor <- cbind(RMA_damage_IPNW_exogenous[2:3], RMA_damage_IPNW_exogenous[5:6], RMA_damage_IPNW_exogenous[8], RMA_damage_IPNW_exogenous[11], RMA_damage_IPNW_exogenous[13:14], RMA_damage_IPNW_exogenous[16:17])
fit <- princomp(RMA_damage_IPNW_exogenous_factor, cor=TRUE)
autoplot(prcomp(RMA_damage_IPNW_exogenous_factor, scale = TRUE, center = TRUE), frame = TRUE, frame.type = 'norm', data = RMA_damage_IPNW_transformed, colour = 'state', loadings = TRUE, loadings.label = TRUE, loadings.label.size = 3, label = FALSE, label.size = 2, legend = FALSE) + theme_bw() + geom_text(vjust=-1, label=rownames(RMA_damage_IPNW_exogenous)) + theme(legend.position = "none") + ggtitle("Principle Components: IPNW insurance loss\nby county with damage cause loadings ") + theme(plot.title = element_text(size =28, family = "Serif")) + theme(axis.title.y = element_text(family = "Serif", size=16), axis.title.x = element_text(family = "Serif", size = 16), axis.text.x = element_text(size=rel(1.5), hjust = 1, family = "Serif"), axis.text.y = element_text(size=rel(1.5), hjust = 1, family = "Serif"))
#pca_RMA_damage_IPNW_exogenous$rotation
#plot(pca_RMA_damage_IPNW_exogenous$rotation)
std_dev <- fit$sdev
pr_var <- std_dev^2
#pr_var[1:10]
prop_varex <- pr_var/sum(pr_var)
plot(prop_varex, xlab = "Principal Component",
ylab = "Proportion of Variance Explained",
type = "b", main = "Scree Plot")
Step 14: PCA: IPNW insurance loss, by YEAR with DAMAGE CAUSE loadings: 2001 to 2015
In Step 5, we perform a PCA for the insurance loss for the inland PNW, for all commodities by year, with damage cause as the factor loadings. Data from 2001 is 2015 is used. We additionally have generated a scree plot that shows the proportion of variance explained by the individual components.
#IPNW damage cause as loadings using years as groups
library(reshape2) ; RMA_damage_IPNW_transformed_year <- dcast(RMA_damage_IPNW, year ~ damagecause, value.var = c("loss"), sum)
#RMA_damage_PNW_transformed_year$statecounty <- paste(RMA_damage_PNW_transformed_year$county, "_", RMA_damage_PNW_transformed_year$state, sep="")
rownames(RMA_damage_IPNW_transformed_year) <- RMA_damage_IPNW_transformed_year$year
#RMA_damage_IPNW_exogenous_year <- RMA_damage_IPNW_transformed_year[,2:28]
RMA_damage_IPNW_exogenous_year <- RMA_damage_IPNW_transformed_year
#IPNW with years as the group
RMA_damage_IPNW_exogenous_year <- cbind(RMA_damage_IPNW_exogenous_year[3:4], RMA_damage_IPNW_exogenous_year[7], RMA_damage_IPNW_exogenous_year[9], RMA_damage_IPNW_exogenous_year[,12:15], RMA_damage_IPNW_exogenous_year[18:19])
fit <- princomp(RMA_damage_IPNW_exogenous_year, cor=TRUE)
autoplot(prcomp(RMA_damage_IPNW_exogenous_year, center = TRUE, scale = TRUE), data = RMA_damage_IPNW_transformed_year, loadings = TRUE, loadings.label = TRUE, loadings.label.size = 3, label = FALSE, label.size = 2, legend = FALSE) + theme_bw() + geom_text(vjust=-1, label=rownames(RMA_damage_IPNW_exogenous_year)) + theme(legend.position = "none") + ggtitle("Principle Components: IPNW insurance loss\nby year with damage cause loadings") + theme(plot.title = element_text(size =28, family = "Serif")) + theme(axis.title.y = element_text(family = "Serif", size=16), axis.title.x = element_text(family = "Serif", size = 16), axis.text.x = element_text(size=rel(1.5), hjust = 1, family = "Serif"), axis.text.y = element_text(size=rel(1.5), hjust = 1, family = "Serif"))
#autoplot(prcomp(RMA_damage_IPNW_exogenous, scale = TRUE, center = TRUE), data = RMA_damage_IPNW_transformed, colour = 'county', loadings = TRUE, loadings.label = TRUE, loadings.label.size = 3, label = FALSE, label.size = 3, legend = FALSE) + geom_text(vjust=-1, label=rownames(RMA_damage_IPNW_exogenous)) + theme(legend.position = "none")#pca_RMA_damage_IPNW_exogenous$rotation
#plot(pca_RMA_damage_IPNW_exogenous$rotation)
std_dev <- fit$sdev
pr_var <- std_dev^2
#pr_var[1:10]
prop_varex <- pr_var/sum(pr_var)
plot(prop_varex, xlab = "Principal Component",
ylab = "Proportion of Variance Explained",
type = "b", main = "Scree Plot")
Step 15: PCA: IPNW insurance loss, by COUNTY with DAMAGE CAUSE loadings: 2015
In Step 5, we perform a PCA for the insurance loss for the inland PNW, for all commodities by year, with damage cause as the factor loadings. Data for just 2015 is used. We additionally have generated a scree plot that shows the proportion of variance explained by the individual components.
#just the IPNW 2015
RMA_damage_IPNW_2015 <- subset(RMA_damage_IPNW, year == 2009)
library(reshape2) ; RMA_damage_IPNW_transformed_2015 <- dcast(RMA_damage_IPNW_2015, county + state ~ damagecause, value.var = c("loss"), sum)
RMA_damage_IPNW_transformed_2015$statecounty <- paste(RMA_damage_IPNW_transformed_2015$county, "_", RMA_damage_IPNW_transformed_2015$state, sep="")
rownames(RMA_damage_IPNW_transformed_2015) <- RMA_damage_IPNW_transformed_2015$statecounty
#RMA_damage_IPNW_exogenous_2015 <- RMA_damage_IPNW_transformed_2015[,3:21]
RMA_damage_IPNW_exogenous_2015 <- RMA_damage_IPNW_transformed_2015
#IPNW 2015
RMA_damage_IPNW_exogenous_2015 <- cbind(RMA_damage_IPNW_exogenous_2015[4:5], RMA_damage_IPNW_exogenous_2015[,8:13], RMA_damage_IPNW_exogenous_2015[,15:16])
fit <- princomp(RMA_damage_IPNW_exogenous_year, cor=TRUE)
#RMA_damage_IPNW_exogenous <- scale(RMA_damage_IPNW_exogenous, scale = TRUE, center = TRUE)
pca_RMA_damage_IPNW_exogenous_2015 <- prcomp(RMA_damage_IPNW_exogenous_2015, scale = TRUE, center = TRUE)
#autoplot(prcomp(RMA_damage_IPNW_exogenous_2015, scale = TRUE, center = TRUE), data = RMA_damage_IPNW_transformed_2015, colour = 'state', loadings = TRUE, loadings.label = TRUE, loadings.label.size = 3, label = FALSE, label.size = 3, legend = FALSE) + geom_text(vjust=-1, label=rownames(RMA_damage_IPNW_exogenous_2015)) + theme(legend.position = "none")
autoplot(prcomp(RMA_damage_IPNW_exogenous_2015, center = TRUE, scale = TRUE), frame = TRUE, frame.type = 'norm', data = RMA_damage_IPNW_transformed_2015, colour = 'state', loadings = TRUE, loadings.label = TRUE, loadings.label.size = 3, label = FALSE, label.size = 2, legend = FALSE) + theme_bw() + geom_text(vjust=-1, label=rownames(RMA_damage_IPNW_exogenous_2015)) + theme(legend.position = "none") + ggtitle("Principle Components: 2015 IPNW insurance loss\nby county with damage cause loadings") + theme(plot.title = element_text(size =28, family = "Serif")) + theme(axis.title.y = element_text(family = "Serif", size=16), axis.title.x = element_text(family = "Serif", size = 16), axis.text.x = element_text(size=rel(1.5), hjust = 1, family = "Serif"), axis.text.y = element_text(size=rel(1.5), hjust = 1, family = "Serif"))
#pca_RMA_damage_IPNW_exogenous$rotation
#plot(pca_RMA_damage_IPNW_exogenous$rotation)
std_dev <- fit$sdev
pr_var <- std_dev^2
#pr_var[1:10]
prop_varex <- pr_var/sum(pr_var)
plot(prop_varex, xlab = "Principal Component",
ylab = "Proportion of Variance Explained",
type = "b", main = "Scree Plot")
Step 16: PCA: IPNW WHEAT loss by YEAR with DAMAGE CAUSE loadings: Temporal Mapping of PC1 and PC2
In Step 16, we summarize our PCA findings for the inland PNW, for wheat by year, using damage cause as the factor loadings. We take PC1 and PC2 loadings and plot those for each year.
library(htmlTable)
#year factanal
#factanal
#fit <- factanal(x = RMA_damage_PNW_exogenous_wheat_year_factor, factors = 3, rotation = "varimax")
#load <- fit$loadings[,1:2]
#plot(load,type="n") # set up plot
#text(load,labels=names(RMA_damage_PNW_exogenous_wheat_year_factor),cex=.7) # add variable names
library(FactoMineR)
result <- FactoMineR::PCA(RMA_damage_PNW_exogenous_wheat_year_factor) # graphs generated automatically
#--create table of PC1 and PC2 loadings
PC_table <- cbind(result$var$coord[,1], result$var$coord[,2])
colnames(PC_table) <- c("PC1", "PC2")
PC_table <- as.data.frame(PC_table)
PC_table$PC1 <- round(PC_table$PC1, 3)
PC_table$PC2 <- round(PC_table$PC2, 3)
htmlTable(PC_table,
cgroup = c("2001-2015"),
n.cgroup = c(2),
caption = "Damage cause variable loadings by year",
align = "c",
rnames = rownames(PC_table),
css.cell = "padding-left: .5em; padding-right: .5em;")| Damage cause variable loadings by year | ||
| 2001-2015 | ||
|---|---|---|
| PC1 | PC2 | |
| Cold Wet Weather | -0.325 | 0.919 |
| Cold Winter | 0.834 | 0.462 |
| Drought | 0.845 | 0.111 |
| Excess Moisture/Precip/Rain | -0.174 | 0.961 |
| Fire | 0.588 | 0.123 |
| Flood | -0.431 | 0.832 |
| Freeze | 0.729 | 0.396 |
| Frost | 0.168 | -0.169 |
| Heat | 0.865 | 0.008 |
| Hot Wind | 0.365 | 0.015 |
#--year end
PC_table_year <- cbind(result$ind$coord[,1], result$ind$coord[,2])
colnames(PC_table_year) <- c("PC1", "PC2")
par(mar=c(6, 6.1, 4, 2.1), family = 'serif', mgp=c(4, 1, 0), las=0)
barplot(t(PC_table_year), beside = T, las = 2, col=c("blue", "red"), cex.axis = 1.5, cex.names = 1.5)
#legend("top", c("PC1","PC2"), pch=15,
#col=c("blue","red"),
#bty="n")Step 17: PCA: IPNW WHEAT loss by COUNTY with DAMAGE CAUSE loadings: Spatial Mapping of PC1 and PC2
In Step 17, we summarize our PCA findings for the inland PNW, for wheat by county, using damage cause as the factor loadings. We take PC1 and PC2 loadings and plot each as a map.
#factanal
#fit <- factanal(x = RMA_damage_IPNW_exogenous_factor, factors = 3, rotation = "varimax")
#load <- fit$loadings[,1:2]
#plot(load,type="n") # set up plot
#text(load,labels=names(RMA_damage_IPNW_exogenous_factor),cex=.7) # add variable names
library(FactoMineR)
result <- FactoMineR::PCA(RMA_damage_IPNW_exogenous_factor) # graphs generated automatically
#--create table of PC1 and PC2 loadings
PC_table <- cbind(result$var$coord[,1], result$var$coord[,2])
colnames(PC_table) <- c("PC1", "PC2")
PC_table <- as.data.frame(PC_table)
PC_table$PC1 <- round(PC_table$PC1, 3)
PC_table$PC2 <- round(PC_table$PC2, 3)
htmlTable(PC_table,
cgroup = c("2001-2015"),
n.cgroup = c(2),
caption = "Damage cause variable loadings by county",
align = "c",
rnames = rownames(PC_table),
css.cell = "padding-left: .5em; padding-right: .5em;")| Damage cause variable loadings by county | ||
| 2001-2015 | ||
|---|---|---|
| PC1 | PC2 | |
| Cold Wet Weather | 0.465 | 0.797 |
| Cold Winter | 0.58 | -0.244 |
| Decline in Price | 0.579 | 0.56 |
| Drought | 0.86 | -0.21 |
| Excess Moisture/Precip/Rain | -0.107 | 0.897 |
| Fire | 0.738 | 0.043 |
| Freeze | 0.913 | -0.176 |
| Frost | 0.945 | -0.083 |
| Heat | 0.929 | 0.079 |
| Hot Wind | 0.882 | -0.163 |
#-map PC loadings
PC1_map <- as.data.frame(result$ind$coord[,1])
PC1_map$NAME <- rownames(PC1_map)
colnames(PC1_map) <- c("PC1", "NAME")
PC1_NAME <- cbind(PC1_map, t(as.data.frame(strsplit(PC1_map$NAME, split='_', fixed=TRUE))))
colnames(PC1_NAME) <- c("PC1", "Statecounty", "NAME", "STATE")
PC1_map <- merge(counties, PC1_NAME, by=c("NAME"), all = FALSE)
PC2_map <- as.data.frame(result$ind$coord[,2])
PC2_map$NAME <- rownames(PC2_map)
colnames(PC2_map) <- c("PC2", "NAME")
PC2_NAME <- cbind(PC2_map, t(as.data.frame(strsplit(PC2_map$NAME, split='_', fixed=TRUE))))
colnames(PC2_NAME) <- c("PC2", "Statecounty", "NAME", "STATE")
PC2_map <- merge(counties, PC2_NAME, by=c("NAME"), all = FALSE)
#PC1 map
#-----
PC1_map$PC1[is.na(PC1_map$PC1)] <- 0
## Make vector of colors for values smaller than 0 (20 colors)
rc1 <- colorRampPalette(colors = c("blue", "lightblue", "white"), space = "Lab")(4)
## Make vector of colors for values larger than 0 (180 colors)
rc2 <- colorRampPalette(colors = c("white", "red"), space = "Lab")(16)
## Combine the two color palettes
rampcols <- c(rc1, rc2)
# pal <- colorNumeric(colorRampPalette(c("blue", "lightblue", "white", "lightpink", "red"))(7), na.color = "#ffffff",domain = eval(parse(text=paste("PC1_map$", "PC1", sep=""))))
pal <- colorNumeric(palette = rampcols, na.color = "#ffffff",domain = eval(parse(text=paste("PC1_map$", "PC1", sep=""))))
map <- leaflet(data = PC1_map, options = leafletOptions(zoomControl = FALSE, zoomSnap = 0.1)) %>%
addProviderTiles(providers$Hydda.Base) %>%
addProviderTiles(providers$Stamen.TonerLines) %>%
addMiniMap(
tiles = providers$Stamen.TonerLite,
position = 'topleft',
width = 100, height = 100,
toggleDisplay = FALSE) %>%
#addControl(title, position = "topleft", className="map-title") %>%
fitBounds(exte[1], exte[3], exte[2], exte[4]) %>% addPolygons(fillColor = ~pal(PC1_map$PC1), fillOpacity = .8, weight = 1, stroke = TRUE, smoothFactor = 0.5, color = "black") %>% addPolygons(data = states, fillOpacity = 0, weight = 4, stroke = TRUE, smoothFactor = 0.5, color = "black") %>%
addLegend(pal = pal, values = na.omit(~PC1), group = "circles", title = paste("PC1 Loadings", sep = ""), na.label = "", position = "topright") %>%
#label = lapply(labs, HTML),
addLabelOnlyMarkers(data = PC1_map, lng = lat_long[,1], lat = lat_long[,2], label = lapply(round(PC1_map$PC1, 2), HTML), labelOptions = labelOptions(style = list("font-family" = "serif"), noHide = TRUE, direction = 'middle', textOnly = TRUE, textsize = "20px", col = "white", offset = c(20, 0), markerOptions(riseOnHover = TRUE)
))
addScaleBar(map, position = c("topright"), options = scaleBarOptions())#PC2 map
PC2_map$PC2[is.na(PC2_map$PC2)] <- 0
## Make vector of colors for values smaller than 0 (20 colors)
rc1 <- colorRampPalette(colors = c("blue", "lightblue", "white"), space = "Lab")(4)
## Make vector of colors for values larger than 0 (180 colors)
rc2 <- colorRampPalette(colors = c("white", "red"), space = "Lab")(16)
## Combine the two color palettes
rampcols <- c(rc1, rc2)
#-----
#pal <- colorNumeric(brewer.pal(11,"RdBu"), na.color = "#ffffff",domain = eval(parse(text=paste("PC2_map$", "PC2", sep=""))))
pal <- colorNumeric(palette = rampcols, na.color = "#ffffff",domain = eval(parse(text=paste("PC2_map$", "PC2", sep=""))))
map <- leaflet(data = PC2_map, options = leafletOptions(zoomControl = FALSE, zoomSnap = 0.1)) %>%
addProviderTiles(providers$Hydda.Base) %>%
addProviderTiles(providers$Stamen.TonerLines) %>%
addMiniMap(
tiles = providers$Stamen.TonerLite,
position = 'topleft',
width = 100, height = 100,
toggleDisplay = FALSE) %>%
#addControl(title, position = "topleft", className="map-title") %>%
fitBounds(exte[1], exte[3], exte[2], exte[4]) %>% addPolygons(fillColor = ~pal(PC2_map$PC2), fillOpacity = .8, weight = 1, stroke = TRUE, smoothFactor = 0.5, color = "black") %>% addPolygons(data = states, fillOpacity = 0, weight = 4, stroke = TRUE, smoothFactor = 0.5, color = "black") %>%
addLegend(pal = pal, values = na.omit(~PC2), group = "circles", title = paste("PC2 Loadings", sep = ""), na.label = "", position = "topright") %>%
addLabelOnlyMarkers(data = PC2_map, lng = lat_long[,1], lat = lat_long[,2], label = lapply(round(PC2_map$PC2, 2), HTML), labelOptions = labelOptions(style = list("font-family" = "serif"), noHide = TRUE, direction = 'middle', textOnly = TRUE, textsize = "20px", col = "white", offset = c(20, 0), markerOptions(riseOnHover = TRUE)
))
addScaleBar(map, position = c("topright"), options = scaleBarOptions())